library(SDMTools)#load necessary libraries
#work.dir = 'C:/Users/alice hughes/Documents/Pakistan_models/2050A1B/overlap/'; setwd(work.dir) #define & set the working directory
work.dir = 'C:/Users/Lauren/Documents/Work backup/SEAbat/new/overlap/';setwd(work.dir)

###01.create a vector of species names
files=list.files(pattern='_overlap.asc') #get a list of files in your working directory that include 'asc'

base.asc=read.asc.gz(files[1]); base.asc[which(is.finite(base.asc))]=0
cur.richness=base.asc
fut.richness=base.asc
stability=base.asc

for (tfile in files) {cat(tfile, '\n')
	overlap=read.asc.gz(tfile);current=overlap;future=overlap
	overlap[which(overlap<3)]=0  #set values less than 3 to zero
	overlap[which(overlap==3)]=1  #set values equal to 3 (where current and future overlap) to 1
	
	current[which(current==2)]=0;current[which(current>0)]=1
	future[which(future==1)]=0;future[which(future>0)]=1
	
	if (is.null(stability)) {stability = overlap;  } #define the output & set everything = 0
	stability = stability + overlap
	if (is.null(cur.richness)) {cur.richness = current;  } #define the output & set everything = 0
	cur.richness = cur.richness + current
	if (is.null(fut.richness)) {fut.richness = future;  } #define the output & set everything = 0
	fut.richness = fut.richness + future
	
	}


zlim=range(c(range(cur.richness,na.rm=T),range(fut.richness,na.rm=T)))#find the limits of current and future richness
prop.cur=stability/cur.richness; #proportion of current richness that remains within current range in the future
#range(prop.cur,na.rm=T) #will give the range
prop.cur[which(prop.cur>1.4)]=1.4;prop.cur[which(prop.cur<0.6)]=0.6 # change the limits to 1.4 and 0.6; must be centred on 1 as 1 = no change
pattern.diff=fut.richness/cur.richness #future richness as a proportion of current
#range(pattern.diff,na.rm=T) #will give the range
pattern.diff[which(pattern.diff>1.4)]=1.4;pattern.diff[which(pattern.diff<0.6)]=0.6 #change the limits
pattern.diff[which(is.nan(pattern.diff))]=1
prop.cur[which(is.nan(prop.cur))]=1

pnts=cbind(x=c(113,116,116,113), y=c(-13,-13,-18.5,-18.5)) #define corner points of legend
cols=colorRampPalette(c('orange','yellow','forestgreen',"#003300"))(100)
stability.cols = colorRampPalette(c('red',"orange","yellow",'forestgreen'))(max(stability,na.rm=T)) #get the colors for the plot
prop.cols=colorRampPalette(c('firebrick',"grey","dodgerblue4"))(100)
out.dir=work.dir
png(paste(out.dir,'future_overlap_richness.png',sep=''), width=dim(cur.richness)[1]*3+30, height=dim(cur.richness)[2]*2+80, units='px', pointsize=20, bg='grey') #start the plot
	par(mfrow=c(2,3),mar=c(2,2,2,2)) #set any plot margins
	image(cur.richness, ann=FALSE,axes=FALSE,col=cols, zlim=zlim) #plot the richness
	legend.gradient(pnts,cols=cols,limits=zlim, title='Richness', cex=1.5)
	mtext('Current',side=1)
	image(fut.richness, ann=FALSE,axes=FALSE,col=cols, zlim=zlim)
	mtext('Future 2085 RCP8.5',side=1)
	plot(10,10,ann=FALSE,axes=FALSE,type='n')#an empty plot
	image(pattern.diff, ann=FALSE,axes=FALSE,col=prop.cols, zlim=c(0.6,1.4))
	legend.gradient(pnts,cols=prop.cols,limits=c(0.6,1.4), title='Proportion', cex=1.5)
	mtext('Proportion of current richness',side=1)
	image(stability, ann=FALSE,axes=FALSE,col=stability.cols, zlim=range(stability,na.rm=T))
	legend.gradient(pnts,cols=stability.cols,limits=range(stability,na.rm=T), title='Richness', cex=1.5)
	mtext('Stability: areas of current and future range overlap',side=1)
	image(prop.cur, ann=FALSE,axes=FALSE,col=prop.cols, zlim=c(0.6,1.4))
	legend.gradient(pnts,cols=prop.cols,limits=c(0.6,1.4), title='Proportion', cex=1.5)
	mtext('Proportion of current richness retained',side=1)
dev.off()
